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FOREWORD 


This  report  is  concerned  with  the  problem  of  real-time  filtering 
of  gravity-gradiometer  readings  on  board  a uniformly  horizontally  moving 
vehicle  to  produce  estimates  of  the  gravity  deflection  and  gravity  anomaly 
at  the  vehicle's  position. 

The  report  is  divided  into  two  chapters.  In  the  first  chapter  the 
gravity  variations  are  modeled  as  arising  from  random  variations  in  mass 
density  along  a line  immediately  below  the  vehicle  at  some  depth.  Using 
first  a discrete  version  of  this  model  in  which  the  state  consists  of 
the  density  at  a moving  sequence  of  position  extending  from  far  behind  to 
far  in  front  of  the  vehicle,  the  covariance  of  the  steady-state  optimal 
filter  is  obtained  by  the  Chandrasekhar  algorithm.  To  confirm  these  re- 
sults, a continuous  model  is  introduced  and  the  optimal  filter  is  obtain- 
able from  Wiener  filter  theory.  The  solutions  of  the  corresponding 
Wiener-Kopf  equations  for  the  filter  transfer  functions  take  relatively 
simple  forms  when  the  measurement  consists  of  the  spatial  derivative,  in 
the  direction  of  motion,  of  the  quantity  to  the  estimated,  and  when  the 
measurement  accuracy  is  sufficiently  high.  The  corresponding  transfer 
functions  in  this  "asymptotic"  limit  are  of  first  or  second  order  depend- 
ing on  whether  the  gravity  anomaly  or  gravity  deflection  is  being  estim- 
ated. 

In  the  second  chapter,  a more  realistic  continuous  model,  due  to 
Heller,  is  used  to  describe  the  gravity  variations.  Asymptotic  Wiener 
filter  results  for  high  measurement  accuracy  are  found  to  take  the  same 
form  as  in  the  first  chapter.  These  results  are  extended  to  cases  in 
which  more  than  one  measurement  is  incorporated  into  an  estimate.  The 
results  concerning  gravity  deflection  estimation  are  unambiguous:  cross- 
track deflection  needs  only  first-order  transfer-functions,  but  in-line 
deflection  needs  second-order  transfer  functions.  The  inclusion  of  the 
gravity-gradient  component  r (x  being  the  direction  of  motion)  in  the 
estimation  of  the  gravity  anomaly  gz  does  not  lead  to  simple  asymptotic 
forms  for  the  filter  transfer  functions.  The  results  obtained  in  Ch.  I 
with  the  discrete  model  suggest,  moreover,  a substantial  increase  in 

-iii- 


1 


accuracy  when  T is  included,  along  with  T , in  the  estimation  of 
■XX  xz 


Chapter  II  concludes  with  the  outline  of  a method,  based  on 
rational  approximation  to  the  transcendental  spectral  density  components, 
for  constructing  a steady-state  filter  without  the  asymptotic  approxima- 
tion used  earlier. 


|K  * V»vr 


CHAPTER  I 


A . INTRODUCTION 


The  disturbance  mass  model  which  causes  the  deflection  of  the 
vertical  and  gravity  anomaly  is  assumed  to  be  a one-dimensional  hori- 
zontal line  mass  distribution  below  the  vehicle's  path.  The  stochastic 
property  of  the  line  mass  density  is  considered  to  be  a white  noise. 

The  intensity  of  the  white  noise  and  the  depth  of  the  line  mass  are 
chosen  to  produce  the  same  root  mean  square  values  and  correlation  dis- 
tance for  the  gravity  deflection  as  the  measured  values  on  the  earth 
surface. 

The  optimal  filter  for  the  gradiometer  measurement  is  sought  by 
two  methods.  First,  we  convert  the  model  into  a discrete  model;  then 
find  the  steady-state  optimal  filter  by  the  Chandrasekhar  algorithm. 

The  other  method  consists  of  finding  the  stationary  filter  by  solving  a 
Wiener-Hopf  equation. 

The  results  obtained  by  both  methods  are  in  good  agreement  with 
each  other.  Of  noteworthy  interest  is  that  the  double  measurement, 
r together  with  - T ) yields  information  on  g comparable 

ZX  XX  ZZ  X 

with  that  from  T alone,  and  on  g much  better  than  that  from  T 

xx  ' xz 


B.  DISCRETE  MODEL 

B-l  System  Equations 

In  a coordinate  frame  fixed  to  the  vehicle  which  moves  with 
velocity  v with  respect  to  the  earth,  the  mass  distribution  of  the 
earth  seems  to  be  time-varying  and  the  rate  of  change  of  the  mass  density 
is  described  by  the  following  partial  differential  equation: 


+ v • Vp  = 0 


(1-1) 


-1- 


where  p denotes  the  mass  density  of  the  earth,  and  V Is  the  gradient 
operator. 

For  simplicity,  we  assume  that  the  disturbance  mass  is  concentrated 
on  a straight  line  with  finite  length  below  the  vehicle  path,  see  Fig. 

1-1,  when  the  vehicle  travels  with  constant  speed  v toward  negative 
x direction,  (l.l)  reduces  to 

- y MSjux1  = o.  ( 1 • 2 ) 

ot  dx 


At  the  boundary  x = -l,  the  density  p(t,  - l)  may  be  considered 
as  a random  process  which  is  assumed  to  be  a white  noise  with  zero  mean 
because  we  are  interested  in  only  deviation  from  the  mean: 


E |p(t,-J)p(t',  -J)|  = | &(t-t') 


(1.3) 


where  E denotes  expectation  operator,  and  q is  the  power  spectral 
density  of  the  white  noise.  Though  (l.2)  does  not  have  any  process  noise, 
the  boundary  condition  (l.3)  always  brings  uncertainty  into  the  system. 

On  the  vehicle,  the  gravity  and  gravity  gradient  due  to  the  dis- 
turbance mass  are  expressed  in  terms  of  integrals  of  the  mass  density 
multiplied  by  weighting  functions: 


where 


= Y j wkp(t,x)dx, 
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(a)  Coordinate  Frame 


(b)  Discrete  Mass  Model 


COORDINATE  FRAME  AND  DISCRETE  MASS  MODEL 


B-2 


Discrete  Formulation 


The  formal  application  of  Kalman  filter  theory  leads  to  a partial 
differential  equation  for  the  states  error  covariance  . To 
avoid  this  difficulty,  here,  the  system  equations  (l.2),  (l.3)  and  ( 1 . 4 ) 
are  discretized  such  that  a distributed  parameter  system  becomes  a 
lumped  parameter  system  and  the  conventional  Kalman  filter  theory  can 
be  applied. 

The  time  and  space  increments,  At  and  Ax,  respectively,  are 
chosen  as  vAt  = Ax  so  that  nondissipative  property  of  the  original 
equation  (l.2)  can  be  preserved).  In  this  case,  (l-2)  may  be  written  as 

= Mj(t^),  j = 1,  N -1  (l.5) 

where  Mj(t^)  denotes  the  mass  of  the  j-th  segment  at  time  t^.  The 
boundary  condition  (1.3)  may  be  replaced  by 

E(Ml(ti)Ml(te)j  = qvAt  6ie  (1,6) 

where  5,  is  the  Kronecker  delta, 
ie 

The  expressions  for  gravity  and  gravity  gradient  (l.4)  may 


easily  be  converted  into  a discrete  form  as  follows: 

N 


B-3  Chandrasekhar  Algorithm 

In  order  to  make  a good  approximation,  we  have  to  take  the  length 
of  the  line  mass,  2 i,  as  large  as  possible,  and  the  space  increment  Ax, 
as  small  as  possible.  When  the  parameters  of  the  system  are  not  time- 
varying  but  constant,  the  Chandrasekhar- type  algorithm  developed  by 
Kailath  et  al.  [Ref.  l]  has  the  possibility  of  substantially  reducing 
the  amount  of  computation  and  compu  r capacity  required  below  those 
necessary  with  the  Riccati  equation.  A brief  explanation  of  the 
Chandrasekhar  algorithm  is  given  in  Appendix  1. 

our  case,  the  state  vector  is  the  mass  of  each  segment 

nr 

[M  (tj)}  , j = 1,  . ..,  N,  and  matrix  F in  App.  1 is  given  by 
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(l.o'l 
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1 0 


The  process  noise  distribution  matrix  G is  given  by 


G — [ 1,  G,  ...,  0]  • 


(1.9) 


The  measurement  matrix  H is  obtained  from  (l.7).  For  example,  when 

T is  used  as  a measurement,  H is  given  by 
xz 


H = <7 


-3Dx  . 

J 


j - 1 , • • • , N 


(1.10) 


The  power  spectral  density  Q of  the  process  noise  is  given  by 


= qvAt  = qAx  . 


(l.ll) 


Q 


The  power  spectral  density  of  the  measurement  noise  R is  given  by 


for  one  measurement  ( 1.12a) 

for  a double  measurement  (l-12b) 

where  r is  the  power  spectral  density  of  the  continuous  measurement 
c 

noise. 

Since  the  dimension  of  process  noise  is  only  one  and  that  of 

measurement  is  one  or  two,  the  number  of  computations  for  one  step  is 
2 2 

on  the  order  of  N x 2 or  N x 3.  Hence,  as  the  number  of  the  segment 
N becomes  large,  the  superiority  of  the  Chandrasekhar  algorithm  over 

3 

the  Riccati  equation  becomes  clearer.  The  latter  needs  an  order  of  N 
computations  for  one  step. 


R = 77 


r 

c 

At 


R = 


c 

At 


— 0 


0 


c 

At 


B-4  Power  Spectral  Density  and  Auto-correlation  Function 

Before  we  proceed  to  numerical  computation  discussed  in  the  pre- 
vious section,  we  mention  the  gravity  field  produced  by  an  infinite  line 
mass  with  white  noise  spectrum.  Extending  the  integral  limits  in  (l.4) 
to  infinity,  (l.4)  may  be  regarded  as  convolution  integrals  of  mass 
density  p and  weighting  functions.  Since  the  Fourier  transform  of 
the  weighting  functions  are  given  by  modified  Bessel  functions  of  the 
second  kind,  the  power  spectral  density  of  the  gravity  and  gravity  gr  id- 
ient  are  easily  obtained.  For  the  deflection  of  the  vertical  along  the 
track,  g^,  and  the  gravity  anomaly  g,  we  find  that 


$ ((‘>)  = 4y2qu)2K2(Du) 
Bx 

$g(w)  = 4y2q<i2K2(Dw) 


( 1 . 13a) 
( 1.13b) 


where  KQ  and  are  zero  and  first-ordt  modified  Bessel  functions, 


f 


1 

; 


respectively.  The  rms  values  of  g^  and  g are  given  by 

(g  ) = £ */nq'2D3  ( 1 . 14a) 

' x rms  2 y ' ' 


(g)rms  = I ^^q/2b3  . (1.14b) 

Auto-correlation  functions  of  g^  and  g are  obtainable  not  analytically 
but  numerically,  and  shown  in  Pig.  1-2. 

Two  parameters  to  be  determined,  namely,  white  noise  intensity  and 

depth  D,  are  chosen  such  that  the  resulting  rms  value  and  correlation 

distance  of  g are  the  sea-level  values  g~  X 8 arcsec,  and  20  n.  mi., 
x 2 u 

respectively  [Ref.  2],  where  gQ  = 9.8  m/s  and  the  correlation  distance 
is  defined  as  the  shift  distance  at  which  the  ACF  drops  to  l/e  of  that 
for  zero  shift.  The  result  is; 

!2  -85,4 

y q = 1.65  X 10  km  /sec 

D = 36  km 

B-5  Numerical  Results 

We  have  investigated  the  estimation  error  covariance  of  the  de- 
flection of  the  vertical  and  the  gravity  anomaly,  when  we  used  different 
components  of  the  gravity  gradient  tensor  as  the  measurements.  The 
results  are  shown  in  Figs.  3 and  4.  The  conclusions  drawn  are  as  follows. 

1.  As  expected,  among  "single"  measurements,  r is  preferable 
to  T for  estimating  its  integral  g while  r is  preferable  to 

ZX  X zx 

r for  estimating  its  integral  g. 

2.  The  double  measurement,  T together  with  %(T  - T ) yields 

ZX  * XX  zz 

information  on  g^  comparable  with  that  from  r alone  and  on  g much 

better  than  that  from  T alone.  These  measurements  consist  of  the  out- 

xz 

puts  of  a single  rotating  gradiometer  with  its  spin  axis  aligned  to  the 
g-axis . 


i 

I 


I 
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FIGURE 


At  this  time,  we  do  not  understand  fully  the  reason  for  Conclu- 
sion 2.  A possible  explanation  is  as  follows. 

For  a one-dimensional  gravity  source,  we  have  an  identity,  namely, 

T = - g/D.  Since  T + P + T is  always  zero,  we  find  that 

yy  ' xx  yy  zz  ’ 

P = — r + g/D.  Therefore,  T has  not  a weak  but  strong  cor- 


zz 


xx 


relation  with  both  g and  g.  Since  — T ) = T - g/2D, 

we  can  say  that,  qualitatively,  the  double  measurement  provides  a 
better  estimate  of  g than  that  from  r alone.  On  the  contrary, 

it  cannot  provide  better  estimate  of  g than  that  from  r 

X XX 

alone.  However,  since  we  have  a good  estimate  of  g,  at  least  for 

relatively  low  measurement  noise,  the  measurement  ~ ^ 

gives  a satisfactory  estimate  of  and  hence  of  g^.  For 

example,  for  1 Eb’tvos  measurement  noise  rms  value  with  10  sec 

average,  the  estimation  error  covariance  of  g/D  is  obtained  to 

be  on  the  order  of  0.001  E.  Assuming  that  the  correlation  time 

is  on  the  order  of  20  mins,  we  can  say  that  the  maximum  power 

-21  -3 

spectral  density  level  is  of  the  order  10  (sec  );  hence, 

much  smaller  than  that  of  the  measurement  error  which  is  of  the 
—17  —3 

order  10  (sec  ). 


The  computer  program  is  shown  in  Appendix  2.  The  computer  language 
used  is  not  FORTRAN  but  APL  (A  Programming  Language). 


B-6 


Wiener  Filter  Theor 


If  the  filtered  estimate  y of  some  stationary  process  y, 
based  on  a noisy  vector  measurement  n^+  n,  where  n is  m-dimensional 
white  noise,  is  given  by: 


y(s)  = if  (s)[z(s)  + n(s)] 


(1.15) 


where  the  transfer  function  vector  f(s)  has  only  left  half  plane  (LHP) 

^ A 

poles,  the  mean-squared  value  of  the  estimation  error,  y = y - y,  is: 


’ 2‘J  I j.  | 


*T(-s)(<J>  (s)  + 4>  )f(s)  - *T( -s) 3>  (s) 

n zy 


*«/’8W8)  + ^yy^  SM  ds  ’ 


where  <Pyy(jw)  is  the  spectral  density  of  y,  the  Fourier  transform 
of  its  ACF  (autocorrelation  f unction) , i .e . , 


/jw 

C (*)  cos  uxdx  , 

y 


(c  (x)  being  the  ACF  of  y(x)),  and  $ , $ are  the  transforms  of  the 

y zz  zy 

ntKm 

ACF  of  the  vector  z and  of  its  cross-correlation  with  y. 

The  optimum  ^(s)  is  determined  by  the  requirements: 


l)  CT~  in  ( 1 . 15)  must  be  finite; 


2)  (*  (J)  + 0*(a)  - 4>  (s)  has  no  LHP  poles. 


The  minimum  mean—  squared  estimation  error  is  then 


<4  = 777  \ - <*>_„( -s)*(»)  3 ds  , (1.18; 


where  c is  any  contour  large  enough 
to  include  all  the  LHP  poles  of  the 
integrand  in  (1.17). 


c 


s-plane 


We  will  apply  the  method  of  the  previous  section,  B-6,  to  the 

estimation  of  g (=  -g)  from  the  single  measurement  T , the  spatial 
Z XZ 

derivative  of  g in  the  direction  of  motion.  The  continuous  analogue 
z 

of  the  second  component  of  equation  (l.7)  is: 


7D 

= 3 © 


(1.19) 


where  P is  linear  density,  a random  process,  and  ^xl  denotes  convo- 
lution.  Now  the  Fourier  transform  of  YD/ r is 


/-<*>  YD  cos  wx  2y 

\ T~2  Xidx  = ~ aKi 

J -<*>  (x  + D V 


(ft)  , 


( 1 .20) 


where  ft  = uD,  and  is  a modified  Bessel  function.  Hence, 

<P_  (Jw)  = ~~2~  ft2K ^(ft)cPp 

BzBz  D 


(1.21) 


and  cpp  and  qjp  r are  obtained  from  (l.2l)  by  multiplication  by 

xz,  qz  xz  xz 

2 

by  -ju>  and  w . <S>  ^ is  the  earlier  q 

To  proceed  further,  we  need  a rational  approximation  to  ftK^ft). 
A suitable  approximation  here  is 


, . . 22.173  

ftK  (ft)  2 2 2 

1 (2.252  + ft  ) 


(1.22) 


The  related  function  fi  k (ft)  and  its  approximant  are  shown  in  Fig.  1-5 
2 1 

(ft  K^(ft)  is  proportional  to  the  transform:  f /P.) 


Introducing  S = SD,  condition  2 of  Eq.  (l,17)  becomes 


V^-s2) 


4B  y s 


4,  2 2.4  + fPn  + 3,  2 2.4 

D (a  - S ) D (a  -S  ) 


(1.23) 


has  no  LHP  poles;  where  B = 22.173,  a = 2.252. 


The  optimal  transfer  function  \|r(s)  is  thus  given  by 


1 - 5 *(s) 


(s  + a)* 

II  < 8 + 

J=i  J 


(1.24) 


where  -Sj(j  = 1>  . ..,  4)  are  the  four  LHP  roots  of  the  equation 


D4Cp  (a2  - S2)4  = 4B2q.  S2  . 

n P 


(1.25) 


Now,  for  small  rp  , Eq.  (l.25)  has  a small  LHP  root: 


S = - 


12  lr p 

D W n 

» m 


(1.26) 


as  well  as  three  large  LHP  roots  given  by 


4B2  q ip 

S6  - — P = 0 . 

D q> 

n 


A low-frequency  approximation  to  (l.24)  is  thus: 


15 
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l - S\|f  (s  ) = 


1 - ♦(«) 


2B7 


D V®n/«fp 


where 


L = 


2TB 


a D ^7^: 


(■*#£) 


_L/L, 

i 

s +I 


(1.27) 


The  result  (1.27)  could  have  been  obtained  more  directly  by  replac- 

4 

ing  0K^(0)  by  B/a  (e?  0.87  instead  of  unity,  which  is  really  the  correct 
limit  as  Q -*  0)  and  requiring  that 


4B 


V 


g 2 S[1  ~ s^(s)]  + *Pn'lf(s ) have  no  LHP  poles 

a D 


1 — s\|r(s)  -»  0 as  s 4 « 


(1.28) 


The  mean-squared  estimation  error  is 


8L 


2*j  £ 


2 2 

4B  ? V 

go  C1  “ st(s)]ds  = 

a D 


4B27V 


8 2 
a D L 


(1.29) 


If,  for  example,  measurement  accuracy  is  1 E every  10  sec  at  speed  50  m/s, 

—18  / 4 2 2 —2 

cp  = 0.5X10  km/sec  , and,  using  D = 36  km,  L = 8700  km,  and  rw  /<T  = 1.1x10  . 

gzgz 

The  straight  line  in  Fig.  1-4  shows  the  asymptotic  behavior  of  this  ratio 
according  to  (1.29).  Note  that  L is  "integration  length"  for  the  first- 
order  integrator"  i|r(s)  = l/(s+  l/L).  The  numerical  value  8700  km  is 
certainly  excessive  for  some  applications. 

We  turn  next  to  the  estimation  of  g^  from  the  single  measurement 

rxx-  Since 


yx  ^ „ 

Bx  = “a  © p - 
r 
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and  75c/ r has  Fourier  transform 


Low  frequency  approximation  to  cp  is 

Kxgx 


4y^-4Bfiy 


V,  mT\Ti9" 


(1.30) 


and  (p  and  cp  are  obtained  by  multiplication  by  -jco  and  oj 

9 ^'■v  * vv  -v-v 


XX'  X 


XX  XX 


The  optimum  filter  transfer-function  is  thus  obtainable  from  the 
requirements 


— 64B27 2qp 


12 


P 3 

s [1  - Si)f(s)]  + q>n\|r(s)  has  no  LHP  poles; 


»[1  - s\|r(s)]  ->  0 as  s -» 


(1.31) 


Hence 


2 

S + L 


\|r(s>  = 


2 2s  2 

■ 


where  (l/4)L  = 


64EfV2qp/; 

12 

a ffi 


(1.32) 


The  mean-squared  estimation  error  is 


8L 


c 

64BVqpn 


64B272s2qp 


12.  3 
a L 


12 


[l  - s^(s)]ds 


(1.33) 


— /Bv"  m3/4  m1/4 
a 


17. 


If,  again,  measurement  accuracy  is  1 Etttvos  every  10  sec  (at  speed 
50  m/sec),  we  find: 


= 720  km 


= 0.26  X 10 


The  straight  line  in  Fig.  1-3  shows  the  asymptotic  behavior  of  this  ratio, 
according  to  (1.33). 
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Chapter  II 


TWO-DIMENSIONAL  GRAVITY  SOURCE  TREATMENT 


A.  THE  HELLER  MODEL 

This  model  assumes  several  independent  layers  below  the  earth  sur- 
face, of  "white-noise"  potential  variations  UWi(i  = 1,  ...  5).  The 
external  gravity  potential  fluctuation  U at  radius  r from  the 
earth's  center  is  expressible  by  Poisson’s  integral  formula 


= E — r-~ f — — ^ 

i 4 n J (r  + R^  - 


U(r , e , cp)  = 


d fit  (1) 


2rR  cos  \|f)‘ 


where  R^  is  the  radius  of  the  i-th  layer,  ,ir  denotes  the  angle  be- 
tween directions  (9,  <p)  and  (0',  cp* ),  and  dfi'  = sin  0'  d0 ' dqp' . 

For  the  description  of  fluctuations  with  wavelengths  substantially 
less  than  earth  radius  a 'flat  earth’  approximation  to  (1)  is  sufficient: 


U(x,y,z) 


■ « 0? 


(z  + Di)Uw  (x',y’) 
i 


[(x-x1)  +( y-y1  ) +(z+D±)  ]‘ 


dx'dy',  (2) 


where  z is  altitude  above  earth  surface,  and  D^  in  the  depth  of  the 
i-th  layer. 

Note  that  formula  (2)  corresponds  not  to  two-dimensional  white-noise 
mass  density  variations  at  depths  D^  but  rather  to  layers  of  random 
mass  density  "doublets": 


±±±  + ± + + ±±  •* 
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Heller  proposes  five  separate  layers  but  only  the  three  shallowest 


layers  contribute  appreciably  to  the  short  wavelength  fluctuations. 

Table  1 gives  the  depths  and  the  spectral  densities  qx 

for  the  potential  variations  U 


Table  1 


i 

D.  (km) 

l 

6,  4. 

<p  (km  /sec  ) 

i 

16.3 

7.13  x 10-8 

2 

92.5 

1.07  X 10-4 

3 

390.5 

1.16  x 10-2 

I I 


The  two-dimensional  Fourier  transform  of  relation  (2)  takes  a 
simple  form 


U[w  , w ; z 
x y 


where 


-Z  . U) 

2 e 1 U (oj  , u ) 

Y w±  * y 


(3) 


z.  = z + D. 

l l 


w 


= Ju?  + w2 
v x y 


and  U r«  . u 1 denotes  the  transform  of  U (x,  y): 
w.L  x yJ 


The  two-dimensional  transforms  of  gravity  fluctuations  g^,  g^,  g^ 


are  obtained  by  multiplying  (3)  by  ju  , j oo  , and  -to  respectively. 

x y 


The  mean-squared  fluctuations  in  g^,  g^,  g^  are  thus  obtainable: 


( V v <.) = i 0 V * 1 1 


-2z.<o  . 

Cp  i dto  dw  . 

e x y 


This  leads  easily  (substituting  to^  = to  cos  a,  to^  = to  sin  a)  to 


CP.  r°o  -2z.t0 
2 1 rr2  r f 1 

g 2 g . 4jt  ] 

y Z 1 o 


3 3 


to  dto  = — Z 


x y 


(2zi)/ 


The  two-dimensional  Fourier  transform  of  T > r , etc  are 

yz  zz 


obtained  from 

(3)  by 

multiplication  by 

-jw  W, 

<oa,  etc. 

dimensional  spectral- 

-density 

matrix  relating  gx 

■ V *Z 

is  thus 

A2 
/ x 

to  to 

X y 

jw  to 

X 

-to  to  to 

x y 

2 

-jto  to 

X 

/ 10  to 

X y 

2 

10 

y 

jto  to 

y 

2 

-to  to 

y 

2 

-JCOytO 

-2z  w 

Z cp  e 1 

i 1 

-Jv 

-j(0  (0 

y 

2 

to 

2 

-Jto  to 

y 

3 

-to 

1 

i 

-WWW 

X y 

2 

-10  (0 

y 

2 

Jto  to 

y 

2 2 
to  to 

y 

3 

Jto  to 

y 

1 

\ jwxw2 

2 

JlO  to 

y 

3 

-U) 

3 

-jto  to 

y 

4 

to 

We  need,  of  course,  the  one-dimenslonal  spectral-density  matrix 
along  the  line  y = 0.  This  is  obtained  by  the  operation 


-21- 


We  thus  obtain,  for  each  layer  i,  and  for  states  8^2^,  g /2zi,  gz/2zi, 


r , r : 

yz  zz 


/h3K1(n) 


-jn K^'(n) 


<p< 


2zi)' 


o n3[Kj*(ft)-K  (n)]  o n4[K^«(fl)-K1(n)] 


jtt3K|(a) 


flV(fi) 


ok'1  » (a) 


\ 


0 n4[K^"(a)-K[(a)]  0 a [kj_«  • » (n)-K* « (n)  ] 


\ 


jfl  K[*(n) 


UK” '(H) 


(s 


where 


( 

^(a)  = j 


-iicosh  u 

e cosh  u du 

0 


a modified  Bessei  function. 

The  gravity-gradient  components  F , F , F , which  are  the 

xx  x y xz 

spatial  derivatives  g , g , g is  the  direction  of  motion,  are 

x y z 

easily  included  by  multiplication  by  _+  Joj i . The  component  F is 

x y y 

obtainable  from  F and  F by  Laplace's  equation, 
xx  zz 

The  zero  correlation  between  g , F , and  the  remaining  states 

y xy 

implies,  of  course,  that  the  optimal  estimate  of  g can  involve 


only  the  measurements  P , P , while  the  optimal  estimates  of  g and 

xy  yz’  6x 

g can  involve  only  the  measurements  P , P , P (or  P ). 
z xx  xz’  zz  yy 

B.  ASYMPTOTIC  RESULTS  FOR  VERY  ACCURATE  MEASUREMENTS 

In  accordance  with  the  results  obtained  with  the  one— dimensional 
density  fluctuation  model,  we  anticipate  that,  in  the  case  of  very 
accurate  measurements,  a low-frequency  approximation  to  the  spectral- 
density  matrix  will  suffice.  The  low-frequency  form  of  (6)  is: 


Hence,  for  several  layers,  the  low-frequency  optical— density  matrices 

for  g , T and  for  g , g , T are,  respectivly: 
y yz  x z zz  ’ J 


(8) 


2k3 

-6k4\ 

-V2 

k2S 

-2k3S 

"6k4 

24kJ 

and 

-k2* 

2S 

~6k4 

/ 

V+21V 

"6k4 

24  k 

D 

where  k = l/it  £ [cp/(2z.)^]  and  s denotes  ju>  . 

J i 1 x 

Example  (1):  Estimation  of  g^  from  (or  of  g^  from  Txz) 


iKs)  (rvv  + n) 
A y 


ey  ~ gy 
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Then,  according  to  Eq.  (1.16)  and  Eq.  (8): 


ST  ^ {«- 


s) (-2k  s + q>  )f(s)  + i|r(-s)(2k  s) 
j n j 


+ (-2kgS)i(r(s)  + 2k3|ds 

j°° 

^ {^gt1  + S'K-sHi-S'Ks)]  + Cpn\(f(-s)i|f(s)}ds 


which  is  minimized  by  choosing  \Jr(s)  so  that 


1 - s\|r(s)  -*  0 as  s 


(<Pn  - 2kgS  )\Jr(s)  has  no  L.H  P poles. 


Hence  i)r(s)  = l/[s+(l/L)],  a "first-order  integrator"  with  "integration 
length"  L = J2 kg/ rp^  . 


The  mean-squared  estimation  error  is 


2 1 


2k  [1  - si|r(s) ]ds  = ^2k  cp  . 

<3  on 


This  result  is  entirely  similar  to  that  obtained  earlier  with  the 
one-dimensional  density  fluctuation  model 

If  we  suppose  again  that  independent  measurements  of  r (or  T ) 

xy  xz 

with  an  accuracy  of  1 Eotvos  are  obtained  every  10  sec  at  a velocity 

— 18  4 

of  50  m/sec,  i.e.,  every  0.5  km,  then  rp^  = 0.5  x 10  km/sec  . If 

the  altitude  z is  zero  , we  find,  with  the  aid  of  Table  1, 

— 6 

<J~  = (T~  = 6. 4 X 10  g,  not  too  different  from  the  <j>~  obtained 

gz  gz  gz 

in  Chap  I.  The  integration  length  L n 7900  km  which  may  again  be  excessive 

for  some  applications!  The  fluctuations  themselves  have  rms  values: 


° = J2  ° = \/3kA  = 4.6  X 10-5g 


(The  three  shallowest  layers  of  the  Heller  model  give  a somewhat 
lower  value  than  assumed  earlier.) 

Example  2:  Estimation  of  gv  from  and  f ( >r  of  g^  from  and 


xy 


r .) 

zz 


If  = ^(sxr^  + nx>  + t2(s)(ry/  + n2)  . 

j«o  f /~2k3s2+rvn 

y -J«  I \ -6k4s 


24k5+qp_ 


+ (~2k3s,  6k4)| 


^(s)' 


^(s)- 


+ 2k„)ds 


1 

2rtj 


j® 

I 

-j 


\ 1 2kg[  l+s^  (-s  Jlfl-s^jCOj]  + tf>n\k1(-s)i)r1(s) 

+ 6k4\)r2(-s)[l  - s,|r1(s)] 

+ 6k4[l+s\|fi(-s)]\|r2(s)  + (24k5+cpnH2(-s)\|r2(s)jds 


This  is  minimized  by  choosing  ^(s),  \|r2(s)  so  that 


1 - s^Cs)  -»  0 as  s _>  o° 


^2<S) 


0 as  s -»  oo 


/ 2 

' -2k  s +cp  6k.  s 

3 ^n  4 


-6k. s 
4 


24k5+cp_ 


has  no  LHP  poles. 
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The  asymptotic  result  for  small  rp^  is: 


♦1(  s> 


s + (-) 


VS) 


_zL 


where  X = ka/4kcL.  and  where  the  integration  length  L is  now: 
4 5 


L = 


^ ^n  \ 4k3k5/ 


The  resulting  mean-squared  estimation  error  is: 


°g  = 2nj  ° f2k3^  “ s+1<s)l  + 6k4Mr2(s)]}d£ 
y c 


•V-4-4)' 


In  our  numerical  example,  the  additional  information  provided  by 

the  additional  measurement  is  not  substantial;  0~  is  reduced  by  7% 

g 

y 

percent , and  L by  15  percent . 


Example  3:  Estimation  of  g^  from  r 


“ R. 


0~ 

g_ 


*(s)(Fxx  + n)  . 


j 00 

= f \|r(-s)(k1s4+cpn)\lr(s)  + k1s31)f(s)-^(-s)k1s3-kls2}ds 

-joo 


1 00 

2tT)  } ( ( — kL s2 ) [ l+s\|r(— s )][1— s>)r(s)]  + qp^xjf (— s ) >)r( s ) } ds  , 


-joo 


which  is  minimized  by  choosing  \|r(s)  so  that 
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i 


s[l  - s\|f(s)]  -»  0 as  s -» 


(klS  + )\|r(s)  has  no  LHP  poles 


Hence, 


\|Ks)  = 


2 

S 4 L 

2 2s  2 

S +~  + 7 


a "second-order  integrator  " with  integration  length  L = ( 4 R ^ / qp^ ) 
The  mean-squared  estimation  error  is 


\ = 25T  ^ <-k1»)[1-«t<s)]ds  =1^1  = (4k1'Pn)l/4  ' 

c L 

This,  again,  is  similar  to  the  result  obtained  with  the  one-dimen- 
sional gravity  fluctuation  model. 

-6 

In  our  numerical  example,  a~  = 3.  X 10  g and  L = 2500  km. 

Kx 


Example  4;  Estimation  of  g^  from  and  fxz  . 


if  = *1(»)(rxx  + V + *2<s)rxz  , 

A ,s4+<p 

-j«  v v 


Vs) 


-2k3s  +rpn/  \ *2(s) 


□ /Vs> 

+(k1s3,  -k2s)( 

\*2<s>, 


= (cont'd  next  page) 


-k^s  fS 


j < -kis  [l+s\|r1(-s)][  l-sf(s)]  + (-8)^(8) 

-J-  v 

“ k2s2't,2(-S^  1-s'*,l(s^  ~ k2s2[  1+s'tf1<~s>]t2(s) 


+ t2(~s)(-2k3s  + <pn)\lr2(s)\ds  , 


which  is  minimized  by  choosing  \|r  (s),  \|r  (s)  so  that 

JL  Z 


s [ 1 - s\|f1(s)]  -»  0 

as 

st2  ( s ) -*  0 

as 

4 

3 v 

/ tjts) 

k s +m 

1 n 

~k2S  \ 

V3 

— 2k  s +cp  / 
3 4n 

\^2<s) 

has  no  LHP  poles 


The  asymptotic  solution,  for  small  rp  , is: 

n 


where 


^(s)  = 


2 

s +T 


2 2s  2 

s 


*2<s> 


2 2s  2 

s + — + “j 
L L 


K / \ 1/4 

and  L = — - ( l — 

*1  \ 2klk3/ 


The  mean— squared  estimation  error  is 
2 1 C r . 2.  , 


= 2Jtj  ^ C-k1s2[l-st1(s)]  - k2s2i)f2(s)}ds 

c 

■?(-^)  (-4)' 


,U'  . t 

c 


X 


Again  the  improvement  due  to  the  additional  information,  in  this 
case  Txz,  is  minor.  In  our  numerical  example  L is  reduced  by  10  per- 
cent and  CL  by  only  5 percent. 

^x 

A similar  analysis  applies  to  the  inclusion  of  a 3rd  measurement, 
r , in  the  estimation  of  g . The  improvement  turns  out  to  be  equally 

ZZ  X 

minor . 

However,  the  inclusion  of  the  measurement,  r . in  the  estimation 

xx 

of  gz  from  r does  not  lead  to  a simple  "slow"  filter.  The  results 

from  the  one-dimensional  density  fluctuation  model  suggest  that  there 

may  be  a substantial  improvement  over  the  first  order  integration  of 

T although  this  may  have  been  related  to  the  identity  V = - b g 
xz'  yy  z 

which  applies  only  to  the  one-dimensional  fluctuation  model. 


C.  RATIONAL  APPROXIMATIONS  AND  OPTIMUM  FILTERING 

The  calculation  of  the  optimal  (Wiener)  filter  without  the  low- 
frequency  approximation  (8)  is  possible  only  after  the  Bessel  functions 
K1(fl)  in  (6)  have  been  approximated  by  rational  functions.  The  simplest 
approximation  for  which  will  serve  our  purpose  is 


K^n)  ~ 


r -i  3 • 

> 1 + i7 

L b^  J 


(9) 


The  function  K (0.)  and  its  approximant  one  shown  in  Fig.  1,  for  the 

* 5 

values  b = 8/3  which  preserves  the  correct  value  for  J (1  K (0)df) 
2 o 1 

and  hence  for  0 

‘ xx 

The  direct  evaluation  of  the  optimum  filter  will  be  laborious, 
especially  in  the  case  of  more  than  one  measurement.  An  alternative 
but  equivalent,  approach  is  described  in  the  next  section. 


FIGURE  II- 


D.  FINITE  STATE  MODEL 


4 

; 

i 


We  will  discuss  how  to  construct  a model  in  which  the  components 
of  g and  of  the  gravity-gradient  tensor  T form  part  of  a firtite 
state  which  obeys  a linear  evolution  in  response  to  a finite  white-noise 
input.  The  optimal  estimates  and  their  accuracy  will  then  be  obtainable 
from  the  steady-state  solution  of  a standard  matrix-Riccati  equation. 

For  each  layer  (i  = 1,  2,  3)  we  build  a model  for  g^  and  T ^ as  follows: 


r 

yz 


FIGURE  I I- 2 


r 

yz 


hJ(S)Wi 


gy/2Zi  = h2(S)Wl  + h3(S)W2 


where  S denotes  JG,  and  W , W are  independent  white-noise  inputs 
with  spectral  density  . The  transfer  functions  h^S) 

must  lead  to  the  rational  approximation,  obtained  with  (9),  to  the 
spectral-density  matrix: 


■ 
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r : /u5[K!"'(ii)-K!' 

yz  / 1 1 

'?zt:  Vfl4  K'”(U)  -K^( 


OO]  SJ  [K^'(fi)  - K 


($})]  S2  [K"(fi)  -K(n)] 


i(«)A 

n / 


This  leads  to: 


h,(s)  = 


P„(-s) 


1 U S\6 


where  P (s)  is  a certain  cubic  with  LHP  zeros; 

«3 


h2(s)  . 


H) 


p2(s2) 


(2  * !)  p3<s> 


where  P2  is  a certain  quadratic  polynomial; 


h3(s  ) , 


P4(±  S) 


(l  * I)  P3(s) 


where  P^  is  a certain  quartic. 


The  model  for  g , g , T is  more  elaborate,  as  shown  in  Fig. 
x7  z7  zz 


. i 


FIGURE  I 1-3 


h4(s)w4 

h5(s)w4  + hfi(s)W5 
h?(s)W4  + hg(s)W5  + hg(s)W(. 


where  again  W4,  W,..,  Wg  are  independent  white-noise  inputs  with  spectral 
density  (Pj/it  (2z^)5,  and  the  transfer-functions  hj(S)(j  =4-9)  lead 
to  the  rational  approximation,  obtained  with  (2.9),  to  the  spectral- 
density  matrix: 


'•  • " — 


This  leads  to: 


(a) 


h4(s) 


*;<-■) 


(>  ♦ !)7 


where  P*(s)  is  a certain  quajrtic  with  LHP  zeros; 


(b) 


Pj(S2)(l  - |) 

“»(8>  = ; — it- 

(l  * j)  PjlS) 


where  P*  is  a certain  cubic; 
o 


(c) 


h6(s)  = 


P6(-S) 


(‘  * i)  -;<*> 


where  P (s)  is  a certain  sextic  with  LHP  zeros; 
6 


(d) 


h?(s)  = 


#,  2W  S\2 

s P2(s  )(1  - 5) 


(l  ♦ «)»:<•> 


where  P is  a certain  quadratic; 


The  total  effect  of  all  three  layers  ( i = 1,  2,  3)  is  obtained 
by  adding  in  parallel  three  appropriate  versions  of  Figures  II-2  or  II-3. 

It  is  planned  to  carry  through  this  finite-state  approximation 

and  the  resulting  optimal  filters,  in  order  to  obtain  an  (approximate) 

check  on  the  range  of  validity  of  the  asymptotic  filters  obtained  for 

high  measurement  accuracy.  It  is  hoped  also  to  investigate  the  asymp- 

totic  form,  if  any,  of  the  estimate  of  g when  f and  T are  both 

Z XX  zx 

measured,  and  to  see  if  the  improvement  obtained  with  T in  the  case 

xx 

of  a one-dimensional  density  fluctuation  model  carries  over  to  the  more 
realistic  model. 


Appendix  1 


CHANDRASEKHAR  ALGORITHM 


According  to  Kailath  [Ref.  l],  the  Chandrasekhar  type  algorithm  is 
described  as  follows. 

Suppose  we  have  a process  with  a known  state-space  model  of  the 


*i  = Hxi+Vi 


x.  , = Fx.  + GU. 

l+l  i i 


i S 0 


where 

Euiui 

= QV 

EViVj  ^ R°ij 

§• 

EU.V* 

= 0 

• 

■j 

Exoxo 

II 

o 

EU.x-  , 0 = EV.x.  . 

Assume  that 

[F,  G,  H,  Q,  R} 

are  constant  matrices  with 

n x n,  n x m,  p X n,  m x m,  p x p. 

Let 

d = fjIqF'  + gqg*  - fttqh' (r+htt0h- - nQ 

and  assume  that  we  can  represent  it  (nonuniquely)  as 


D - Wo 

where  LQ  and  MQ  are  n x a and  a X a matrices,  a = rank  D,  and 

Mq  = diag  [1,  1 1,  -1,  -1 -1} 

with  as  many  + l's  as  D has  _+  eigenvalues.  Then  the  quantities 
[K. , R€}  appearing  in  the  estimator  formula 


1 


‘1+1  1 


;1|1-1  * Kl‘"!>'l(>’l  - «i?l|l-l> 


can  be  computed  via  the  equations 


R*  - HLt  (r[)_1L^H' 


= R^  - IjJH*  (Rp^HL. 
= Kt-  FL1(r^)_1L^H' 


= FL.-  Ki(R^)"1HLi 


where  the  initial  values  are 


Ro  = R + hiioh-,  k0=fh0h-,  *rQ  = -m'1. 


The  number  of  computations  required  to  go  from  index  i to  index 

2 *3 

i+1  can  be  seen  to  be  ©(n  (p+Ct))  as  compared  to  if  the 

Riccati  equation  is  used.  For  special  initial  conditions,  matrices 
D,  Lq,  and  are  given  as  follows: 

(a)  HQ  = 0 (perfect  a priori  knowledge  of  the  states); 

D = GQG',  Lq  = GQ^,  and  MQ  = I . 

(b)  IIq  = FI  (stationary  process) 

f“IF'  + GQG'  = K 
D = — TtH*  ( R + HltH’ )-1H~I 

lq  = Tth'(r  + hhh*  )_T/2 

M0  = -!  . 


If  the  error  covariance  matrix  Is  desired,  it  may  be  com- 


puted as 


-38- 
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